//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   testMoment_to_Discrete.cc
 * \author Jeremy Roberts
 * \date   Jul 13, 2011
 * \brief  
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#include "test/Testing.hh"

#include "tran/Traits.hh"

#include "angle/AngleTraits.hh"
#include "angle/Moment_to_Discrete.hh"
#include "angle/Spherical_Harmonics.hh"
#include "angle/Definitions.hh"
#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"
#include "../utilities/SP.hh"

using namespace std;
using namespace slabtran;
using namespace util;

int test1d_1();
int test2d_1();
int test3d_1();

int main(int argc, char *argv[])
{
  Require(argc==2);
  int test = atoi(argv[1]);
  cout << " ------ testMoment_to_Discrete ------ " << endl;
  int testvalue = 0;
  int number_of_tests = 4;
  if (test == 1 || test == number_of_tests) testvalue += test1d_1();
  if (test == 2 || test == number_of_tests) testvalue += test2d_1();
  if (test == 3 || test == number_of_tests) testvalue += test3d_1();
  if (test > number_of_tests || test < 1)
  {
    cout << " test " << test << " not available. " << endl;
  }
  if (testvalue)
    cout << " failed with value = " << testvalue << endl;
  else
    cout << " passed." << endl;
  return testvalue;
}

//===========================================================================//
// 1-D TESTS
//===========================================================================//
int test1d_1()
{
  // Setup the angular system.  Manually set up quadrature for
  // testing to ensure we understand the underlying organization.
  int legendre_order = 3;   // max Legendre order implemented
  int sn_order = 4;         // i.e. 4 directions total
  int number_per_oct = 2;   //
  int number_oct = 2;       // 1-d = left and right

  // Reference angles (first octant)
  double mu[] =
  { 0.8611363115940525752239465, 0.3399810435848562648026658 };

  // Make the test angular mesh.
  AngularMesh<_1D>::SP_angularmesh angular_mesh;
  angular_mesh = new AngularMesh<_1D> (sn_order, GAUSSLEG, false);

  // Make moment-to-discrete operator.
  Moment_to_Discrete<_1D> M(legendre_order);
  M.build(angular_mesh);

  // Designators for the two octants in 1-d.
  int oct[] =
  { 1, -1 };

  for (int o = 0; o < 2; o++)   // octant
  {
    for (int n = 0; n < 2; n++) // angles in octant
    {
      int angle = n + o * number_per_oct;
      for (int l = 0; l <= legendre_order; l++) // moments
      {
        int moment = Moments<_1D>::index(l);
        double factor = 0.5 * (2.0 * l + 1.0);
        double ref = factor * Spherical_Harmonics::Y_lm(l, mu[n]*oct[o]);
        // Test (angle, moment) interface
        TEST(soft_equiv(M(angle, moment), ref));
        // Test (angle, l, m) interface
        TEST(soft_equiv(M(angle, l, 0), ref));
        // Test (octant, angle_octant, l, m) interface
        TEST(soft_equiv(M(o, n, l, 0), ref));
      }
    }
  }
  cout << " ... passed 1-d moment-to-discrete test. " << endl;
  return 0;
}

//===========================================================================//
// 2-D TESTS
//===========================================================================//
int test2d_1()
{
  // Setup the angular system.  Manually set up quadrature for
  // testing to ensure we understand the underlying organization.
  int legendre_order = 3;               // max Legendre order implemented
  int number_moments = (3+1)*(3+2)/2;   // (N+1)(N+2)/2 for 2-d
  int sn_order = 4;                     // 4(6)/8 = 3 per octant.
  int number_per_oct = 3;               //
  int number_oct = 4;                   // 2d = 4 octants.

  // reference angles (first octant)
  double mu[] =
  { 0.350021174581540677777041, 0.35002117458154067777704,
    0.868890300722201205229788 };
  double eta[] =
  { 0.350021174581540677777041, 0.868890300722201205229788,
    0.350021174581540677777041 };
  double xi[] =
  { 0.868890300722201205229788, 0.350021174581540677777041,
    0.350021174581540677777041 };

  // Make the test angular mesh.
  AngularMesh<_2D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_2D> (sn_order, LEVELSYM, false);

  // Make moment-to-discrete operator.
  Moment_to_Discrete<_2D> M(legendre_order);
  M.build(angularmesh);

  // Designators for the four octants in 2-d.
  double oct[4][3] = {{1,1,1}, {-1,1,1}, {1,-1,1}, {-1,-1,1}};

  for (int o = 0; o < number_oct; o++)
  {
    for (int n = 0; n < number_per_oct; n++)
    {
      int angle = n + number_per_oct * o;
      for (int l = 0; l <= legendre_order; l++)
      {
        double factor = (2.0 * l + 1.0) * inv_four_pi;
        for (int m = -l; m <= l; m += 2)
        {
          cout << " l=" << l << " m=" << m << endl;
          int moment = Moments<_2D>::index(l, m);
          double ref = factor * Spherical_Harmonics::
            Y_lm(l, m, mu[n]*oct[o][0], eta[n]*oct[o][1], xi[n]*oct[o][2]);
          // Test (angle, moment) interface
          TEST(soft_equiv(M(angle, moment), ref));
          // Test (angle, l, m) interface
          TEST(soft_equiv(M(angle, l, m), ref));
          // Test (octant, angle_octant, l, m) interface
          TEST(soft_equiv(M(o, n, l, m), ref));
        }
      }
    }
  }

  cout << " ... passed 2-d moment-to-discrete test. " << endl;
  return 0;
}

//===========================================================================//
// 3-D TESTS
//===========================================================================//
int test3d_1()
{
  // Setup the angular system.  Manually set up quadrature for
  // testing to ensure we understand the underlying organization.
  int legendre_order = 3;               // max Legendre order implemented
  int number_moments = (3+1)*(3+1);     // (N+1)^2 for 3-d
  int sn_order = 4;                     // 4(6)/8 = 3 per octant.
  int number_per_oct = 3;               //
  int number_oct = 8;                   // 3d = 8 octants.

  // Reference angles (first octant)
  double mu[] =
  { 0.350021174581540677777041, 0.35002117458154067777704,
    0.868890300722201205229788 };
  double eta[] =
  { 0.350021174581540677777041, 0.868890300722201205229788,
    0.350021174581540677777041 };
  double xi[] =
  { 0.868890300722201205229788, 0.350021174581540677777041,
    0.350021174581540677777041 };

  // Make the test angular mesh.
  AngularMesh<_3D>::SP_angularmesh angular_mesh;
  angular_mesh = new AngularMesh<_3D> (sn_order, LEVELSYM, false);

  // Make moment-to-discrete operator.
  Moment_to_Discrete<_3D> M(legendre_order);
  M.build(angular_mesh);

  // Designators for the eight octants in 3-d.
  double oct[8][3] = {{ 1, 1, 1}, {-1, 1, 1}, { 1,-1, 1}, {-1,-1, 1},
                      { 1, 1,-1}, {-1, 1,-1}, { 1,-1,-1}, {-1,-1,-1}};

  for (int o = 0; o < number_oct; o++)
  {
    for (int n = 0; n < number_per_oct; n++)
    {
      int angle = n + number_per_oct * o;
      for (int l = 0; l <= legendre_order; l++)
      {
        double factor = (2.0 * l + 1.0) * inv_four_pi;
        for (int m = -l; m <= l; m++)
        {
          int moment = Moments<_3D>::index(l, m);
          double ref = factor * Spherical_Harmonics::
            Y_lm(l, m, mu[n]*oct[o][0], eta[n]*oct[o][1], xi[n]*oct[o][2]);
          // Test (angle, moment) interface
          TEST(soft_equiv(M(angle, moment), ref));
          // Test (angle, l, m) interface
          TEST(soft_equiv(M(angle, l, m), ref));
          // Test (octant, angle_octant, l, m) interface
          TEST(soft_equiv(M(o, n, l, m), ref));
        }
      }
    }
  }

  cout << " ... passed 3-d moment-to-discrete test. " << endl;
  return 0;
}

//---------------------------------------------------------------------------//
//              end of testMoment_to_Discrete.cc
//---------------------------------------------------------------------------//
